//========================================================================================
// MerHash.cc
// ----------------------------------------------------------------------------
// $Id: //CBT-depot/CBT++/src/MerHash.cc#1 $
// ------------------------------------------------------------------------------
// $Log: $
//
//========================================================================================
//.......1.........2.........3.........4.........5.........6.........7.........8.........9
//

#include "MerHash.hh"


//--------------------------------------
// * ShiftIn
//--------------------------------------
// Don't call bases twice in same function call because it uses 
// a static char array.
// ShiftIn(nmer, reverseMer, Letter, ngood, proc);
void 
MerHash::ShiftIn(Word &w, Word &r,Word b,int &g,int &p) 
{
  w <<= (BITSPERBASE + BITSUNUSED);  
  w >>= BITSUNUSED; 
  w |= b;

  // KJD FIX This doesn't work with my encoding. Must fix all 
  // Reverse complement computations before switching to my encoding. 
  r >>= BITSPERBASE; 
  r |= (~b & BASEMASK) << (BITSPERWORD - (BITSPERBASE + BITSUNUSED)); 
  g++; 
  
  // Have we gotten a word's worth of bits yet?  If 
  // so, add the smaller of word/revcomp(word) to the 
  // hash.   
  if (g >= (BASESPERWORD-BASESUNUSED)) { 
    Word tw; 
    tw = (w <= r)? w : r; 
    AddCount(tw, p); 
  } 
}

//--------------------------------------
// * Add
//--------------------------------------
// Adds mers from standard input stream. 
void
MerHash::Add(istream &fin)
{
  printf("Hash size is: %ld\n", mHashSize);
  Word nmer = 0;
  Word rev = 0;
  char pc, c;
  int ngood = 0;
  State s = FASTA;
  while ((c = fin.get()) != EOF) {
    if (s == DESCRIP) { 
      if (c == '\n') {
	s = FASTA;
	ngood = 0;
      }
    }else { // s == FASTA
      switch (c) {
	case '>': {
	  s = DESCRIP;
	  mReads++;
	  // A little something for the folks at home. 
	  if (! (mReads & 0xFFFF)) {
	    cerr<<"\r "<<mReads<<" reads\t"<<mNmers<<" nmers\t";
	    cerr<<mDistinct<<"distinct in hash ("<<((mDistinct*100UL)/mNmers)<<"%)"<<endl;
	  }
	  break;
	}
	case 'a': case 'A': {
	  ShiftIn(nmer, rev, A, ngood, mProc);
	  break;
	}
	case 'c': case 'C': {
	  ShiftIn(nmer, rev, C, ngood, mProc);
	  break;
	}
	case 'g': case 'G': {
	  ShiftIn(nmer, rev, G, ngood, mProc);
	  break;
	}
	case 't': case 'T': {
	  ShiftIn(nmer, rev, T, ngood, mProc);
	  break;
	}
	case ' ': case '\t': case '\n': case '*': {
	  // Ignore
	  break;
	}
	default: {
	  nmer = 0;
	  rev = 0;
	  ngood = 0;
	}
      }
    }
  }
}


//--------------------------------------
// * CString2Mer
//--------------------------------------
// 
Word
MerHash::CString2Mer(const char * MerStr)
{
  Word nmer =0;
  int slen = strlen(MerStr);
  for(int i = 0;i< slen;i++){
    char c = MerStr[i];
    switch (c) {
      case 'a': case 'A': {
	ShiftInMer(nmer, A);
	break;
      }
      case 'c': case 'C': {
	ShiftInMer(nmer, C);
	break;
      }
      case 'g': case 'G': {
	ShiftInMer(nmer, G);
	break;
      }
      case 't': case 'T': {
	ShiftInMer(nmer,T);
	break;
      }
      case ' ': case '\t': case '\n': case '*': {
	// Ignore
	break;
      }
      default: {
	nmer = 0;
      }
    }
  }
  return(nmer);
}


//--------------------------------------
// * String2Mer
//--------------------------------------
// 
Word
MerHash::String2Mer(string &MerStr)
{
  Word nmer;
  int slen = MerStr.size();
  for(int i = 0;i< slen;i++){
    char c = MerStr[i];
    switch (c) {
      case 'a': case 'A': {
	ShiftInMer(nmer, A);
	break;
      }
      case 'c': case 'C': {
	ShiftInMer(nmer, C);
	break;
      }
      case 'g': case 'G': {
	ShiftInMer(nmer, G);
	break;
      }
      case 't': case 'T': {
	ShiftInMer(nmer,T);
	break;
      }
      case ' ': case '\t': case '\n': case '*': {
	// Ignore
	break;
      }
      default: {
	nmer = 0;
      }
    }
  }
  return(nmer);
}

//--------------------------------------
// * Mer2CString
//--------------------------------------
// 
char *
MerHash::Mer2CString(Word w) {
  //  static char b[BASESPERWORD+1];
  //	char *b = new char(BASESPERWORD+1);

  char base = '\0';
  int basebits;
  int i;
  unsigned long long bugNum = (BASESPERWORD - BASESUNUSED);

  // This looks buggy because the last value assigned to "base" is not
  // used.  That is an extra assignment with no effect; all the bases
  // and null are properly stored in the (1 + BASESPERWORD - BASESUNUSED) iterations.
  for (i = (BASESPERWORD - BASESUNUSED); i >= 0; i--) {
    tempbuf[i] = base;
    basebits = w & BASEMASK;
    w >>= BITSPERBASE;
    switch (basebits) {
      case A: 
	base = 'a';
	break;
      case C: 
	base = 'c';
	break;
      case G: 
	base = 'g';
	break;
      case T: 
	base = 't';
	break;
    }
  }
  return tempbuf;
}

//--------------------------------------
// * AddCount
//--------------------------------------
//
void 
MerHash::AddCount(Word w, int proc) {
  Index probe, i;
  Count t;
  Word temp;
  Count bitcount = 0;
  
  //#ifndef UNIBATCH
  //for (temp = lo ^ up; temp; temp >>= 5) {
  //	 bitcount ^= temp & 0x1F;
  // }
  //if (bitcount != proc) {
  //	 return;
  // }
  //#endif

  mNmers++;
  
  findslot(probe,w);
  // fprintf(stderr, "%s %ld\n", Mer2CString(w), probe);
  if (! mCount[probe]) {
    mDistinct++;
    mCount[probe] = 1;
    mHash[probe] = w;
  }else if (mHash[probe] != w) {
    fprintf(stderr, "\nFull hash on word %s, probe %ld, insertion %ld\n", 
	    Mer2CString(w), probe, mDistinct); 
    exit(1);
  }else if (((Count)(t = mCount[probe] + 1)) > 0) {
    mCount[probe] = t;
  }
}

//--------------------------------------
// * ShiftInMer
//--------------------------------------
//
void 
MerHash::ShiftInMer(Word &w, Word b) 
{
  w <<= (BITSPERBASE + BITSUNUSED);  
  w >>= BITSUNUSED; 
  w |= b;
  //cout<<"ShiftInMer :"<<w<<endl;
}

//char basebits = (MerStr[i] & 0x06) >> 1;
//shiftinmer(w,basebits);


//--------------------------------------
// * Add
//--------------------------------------
// 
void 
MerHash::Add(string &MerStr,int count){
  //Word w = String2Mer(MerStr);  // Replace with inline to try performance. 

  Word w;
  int slen = MerStr.size();
  for(int i = 0;i< slen;i++){
    char c = MerStr[i];
    switch (c) {
      case 'a': case 'A': {
	shiftinmer(w, A);
	break;
      }
      case 'c': case 'C': {
	shiftinmer(w, C);
	break;
      }
      case 'g': case 'G': {
	shiftinmer(w, G);
	break;
      }
      case 't': case 'T': {
	shiftinmer(w,T);
	break;
      }
      default: {
	cerr<<"ILLEGAL CHARACTER: "<<c<<" in mer attempted add to hash."<<endl; 
	return;
      }
    }
  }
  Add(w,count);
}

//--------------------------------------
// * Add
//--------------------------------------
// 
void 
MerHash::Add(char* MerStr,int count){
  //Word w = String2Mer(MerStr);  // Replace with inline to try performance. 

  Word w;
  //int slen = MerStr.size();
  int stop = 32 - BASESUNUSED;
  for(int i = 0; i < stop; i++){
    char c = MerStr[i];
    switch (c) {
      case 'a': case 'A': {
	shiftinmer(w, A);
	break;
      }
      case 'c': case 'C': {
	shiftinmer(w, C);
	break;
      }
      case 'g': case 'G': {
	shiftinmer(w, G);
	break;
      }
      case 't': case 'T': {
	shiftinmer(w,T);
	break;
      }
      default: {
	cerr<<"ILLEGAL CHARACTER: "<<c<<" in mer attempted add to hash."<<endl; 
	return;
      }
    }
  }
  // cout<<"Adding " << MerStr << "; w = " << Mer2CString(w) << ", count = " << count <<endl;
  Add(w,count);
}



//--------------------------------------
// * Add
//--------------------------------------
// 
void 
MerHash::Add(Word w,int count){

  mTotalInsertions++;

  Index probe,i;
  Count t;
  Word temp;
  Count bitcount = 0;
  const Count maxCount = (Count) 0xFFFF;
  
  findslot(probe,w);

  // fprintf(stderr, "%s %ld\n", Mer2CString(w), probe);

  if (mCount[probe] && mHash[probe] != w) {
    fprintf(stderr, "\nFull hash on word %s, probe %ld, insertion %ld\n", 
	    Mer2CString(w), probe, mDistinct); 
    exit(1);
  }
  else {
    // Not already in the hash. or in hash with different count
    mHash[probe] = w;
    mCount[probe] = (count > maxCount)? maxCount : count;
  }
}

//--------------------------------------
// * Mer2Word
//--------------------------------------
// Mer is in the wrong encoding, so need to translate it to 
// a Word and call Contains on that.  This is just a stop gap 
// way to do this, in order to get something working quickly. 
// Efficiency will be addressed by changing hash encoding. 
//
// KJD Shoot me.  I can't believe how I encoded this. 
// acgtacgtacgt ggggccccttttaaaa
// 0123....
// 
// So given an input sequence or string, the numbering will
// go from left to right as above.  The leftmost base will be zero. 
// My mer encoding puts bases into the upper or lower 32 bits
// depending on whether the base index i < 16 or not.  If < 16, 
// it goes into the Lower32, if > 16, it goes into the upper. 
// Worse, I shift the bases in from right to left.  The result, 
// the bits of the Mer correspond to:
// 
// aaaattttccccgggg tgcatgcatgcatgca
// 
// Horrid.   I don't know what I was thinking, but for now
// it will be easier to have this one function here handle the
// weirdness than it will be to change all of the places in 
// MerHashMap where the encoding is crazy.   At least I think it
// will be....
//
Word
MerHash::Mer2Word(Mer M)
{
  const unsigned char transtable = // 0xb4; // 10110100 (11 -> 10, 10 -> 11, etc)
    (A << (Akjd * 2)) | (C << (Ckjd * 2)) | (G << (Gkjd * 2)) | (T << (Tkjd * 2));
  // printf("%x\n", transtable);
  Index upper32,lower32;
  upper32 = M.mUpper32;
  lower32 = M.mLower32;

  // cout << "Upper32: " << hex << upper32 << endl;
  // cout << "Lower32: " << lower32 << dec << endl;
  Word w = 0;

  Index phhup = 0;
  Index phhlo = 0;
  Index revup = 0;
  Index revlo = 0;
  char kjdbits, phhbits;

  // OK, this is crazy.  My Lower32 bits will become
  // the upper 32 bits for the hash, and as I shift my 
  // bits out to the right, they are added to the left. 

  // Shift in the lower 32 bits (2 bits at a time). 
  // Set up size for divisions
#define BITSPERINDEX (sizeof(Index) * 8)

  for(int i = 0;i< 16;i++){
    kjdbits = (lower32 & 0x00000003) << 1; // grab the lower two bits, * 2
    lower32 = lower32>>2; // Shift them out. 

    phhbits = (transtable >> kjdbits) & 0x3;
    phhup <<= BITSPERBASE;
    phhup |= phhbits;

    revlo >>= BITSPERBASE;
    revlo |= ((Index) phhbits) << (BITSPERINDEX - BITSPERBASE);

    // Shift in the upper 32 bits (2 bits at a time). 

    kjdbits = (upper32 & 0x00000003) << 1; // grab the lower two bits, * 2. 
    upper32 = upper32>>2; // Shift them out. 
	 
    phhbits = (transtable >> kjdbits) & 0x3;
    phhlo <<= BITSPERBASE;
    phhlo |= phhbits;

    revup >>= BITSPERBASE;
    revup |= ((Index) phhbits) << (BITSPERINDEX - BITSPERBASE);
  }
  w = (((Word) phhup) << 32) | phhlo;
  Word r = ~((((Word) revup) << 32) | revlo);
  if (BASESUNUSED) {
    w >>= BITSUNUSED;
    r <<= BITSUNUSED;
    r >>= BITSUNUSED;
  }
  // cout<<"w= "<< hex << w << " " << Mer2CString(w) <<endl;
  // cout<<"r= "<< r << " " << Mer2CString(r) <<endl;
  w = (w <= r)? w : r;
  // cout<<"min= "<< w << dec << " " << Mer2CString(w) <<endl;
  return(w);
}

//--------------------------------------
// * Contains
//--------------------------------------
// More efficient would be to build M and rcM at the 
// same time, and pass both in.  
// 
// 
bool
MerHash::Contains(Mer M)
{
  Word w = Mer2Word(M);	// now returns canonical Mer (lesser of fwd, rev)
  if (ContainsFwd(w)) return(true);
  else return(false);
  //    else{
  //  	 w = RevComp(w);
  //  	 if (ContainsFwd(w)) return(true);
  //  	 else return(false);
  //    }
} 

//--------------------------------------
// * Contains
//--------------------------------------
//
bool
MerHash::Contains(Mer M,int &count){
  Word w = Mer2Word(M);	// now returns canonical Mer (lesser of fwd, rev)
  Count ccount = 0;
  if (ContainsFwd(w,ccount)){
    count = (int) ccount;
    return(true);
  }
  else {
    count = 0;
    return(false);
  }
  //    }else{
  //  	 w = RevComp(w);
  //  	 if (ContainsFwd(w,ccount)){
  //  		count = (int) ccount;
  //  		return(true);
  //  	 }else{
  //  		count = 0;
  //  		return(false);
  //  	 }
  //    }
}

//--------------------------------------
// * ContainsFwd
//--------------------------------------
//
bool
MerHash::ContainsFwd(const char *MerStr)
{
  Word Mer = CString2Mer(MerStr);
  return(ContainsFwd(Mer));
}



//--------------------------------------
// * ContainsFwd
//--------------------------------------
//
bool
MerHash::ContainsFwd(Word w)
{
  Index probe,i;
  Count t;
  Word temp;
  Count bitcount = 0;
      
  findslot(probe,w);

  if ((mHash[probe] == w) && (mCount[probe] > 0)) return(true);
  else return(false);  
}


//--------------------------------------
// * Contains
//--------------------------------------
//
bool
MerHash::ContainsFwd(Word w,Count &count)
{
  Index probe,i;
  Count t;
  Word temp;
  Count bitcount = 0;
  
  findslot(probe,w);

  if ((mHash[probe] == w) && (mCount[probe] > 0)){
    count = mCount[probe];
    return(true);
  }else {
    count = 0;
    return(false); 
  } 
}
//--------------------------------------
// * RevComp
//--------------------------------------
//
Word 
MerHash::RevComp(Word w) {
  Index i;
  Word rev = 0;
  
  // The word w is treated as a string of 2-bit bases
  // (A == 01, T == 11, C == 01, G == 10)
  // So long as we reverse the string basewise,
  // we just need to bit-complement to get the final answer.  
  for (i = 1; i <= BASESPERWORD; i++) {
    //  Shift result left and add in last base of input; then advance input
    rev = (rev << BITSPERBASE) + (w & BASEMASK);
    w >>= BITSPERBASE;
  }
  return ~rev; // bitwise complement == basewise complement 
}

void 
MerHash::SetMerSize(int merSize){
  BASESUNUSED = 32 - merSize; 
  BITSUNUSED=(BASESUNUSED*BITSPERBASE);
}
